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BLOCK  19:  ABSTRACT 

A  method  is  presented  for  computing  three-dimensional  seismic  wave  scattering  from  a 
rough  interface.  It  is  derived  using  a  perturbation  approach  which  requires  interface  height 
perturbations  to  be  small  relative  to  the  wavelengths  of  scattered  waves,  and  interface  slope 
perturbations  to  be  much  less  than  one.  These  validity  conditions  are  based  on  an  order 
of  error  analysis  of  the  truncation  of  the  perturbation  series.  Comparison  of  reflection  and 
transmission  coefficients  generated  by  the  perturbation  method  with  those  generated  from  a 
second-order  finite  difference  method  for  several  rough  interface  models  with  Gaussian  auto¬ 
correlation  functions  shows  that  the  perturbation  method  is  accurate  for  RMS  interface* 
height  deviations  of  less  than  about  a  tenth  of  the  smallest  wavelength  in  the  scattered  field. 
This  result  is  independent  of  RMS  interface  slope  in  the  tested  range  of  0.037  to  0.99. 

Three-dimensional  scattering  results  are  presented  for  normally  incident,  planar  P  and 
SV  waves  on  a  rough  interface  with  a  Gaussian  auto-correlation  function  (RMS  height  is 
0.125  Sj  wavelengths  and  RMS  slope  in  0.03).  For  the  SV  source,  the  scattering  coefficients 
for  the  generated  P,  SV,  and  SH  waves  had  maxima  in  the  azimuth  where  the  scattered  wave 
particle  motion  coincided  with  the  source  particle  motion.  The  maximum  of  the  SV  reflection 
scattering  coefficient  is  37  percent  of  the  mean  planar  interface  reflection  coefficient.  The  P 
source  generated  P  and  SV  waves  with  azimuthally  symmetric  scattering  coefficients.  The 
maximum  of  the  P  scattering  coefficient  in  this  case  is  27  percent  of  mean  planar  interface 
reflection  coefficient. 
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Preface 

The  following  document  is  a  portion  ot  the  Pii.D.  Thesis  of  Michael  P range  entitled 
‘Perturbation  Approximation  of  3-D  Seismic  Scattering”. 


INTRODUCTION 


The  presence  of  a  rough  interface  can  strongly  affect  seismic  waxes  reflected  from  and  trails 
initted  through  that  interface,  even  when  the  scale  of  roughness  is  much  less  than  a  wave 
length.  These  effects  include  changes  in  the  amplitude,  scattering  angle,  frequency  content . 
and  wave-type  conversion  of  the  scattered  wave.  Available  exact  solutions  take  the  form  c! 
integral  eouations  (DeSanto  and  Brown,  1986)  or  finite  difference/finite  element  formula) ions 
(Levander  and  Hill,  1985)  that  are  prohibitively  expensive  to  solve  in  three  dimensions. 

In  this  paper.  I  present  a  perturbation  approach  to  the  solution  of  the  !  hree-dimeusi<  •u  d 
elastic  wave  equation  which  satisfies  welded  boundary  conditions  at  a  rough  interlace.  !  re¬ 
solution  is  an  extension  of  two-dimensional  surface  wave  scattering  formula!  iuusf  Kenne! : . 
1972;  Gilbert  and  Knopoff.  1960).  The  method  height  and  slope  of  interface  irregularities  m 
be  small  with  respect  to  the  wavelengths  of  the  elastic  waves  present.  In  order  to  understand 
the  effects  of  rough  interface  scattering  on  body  waxes,  we  then  present  a  method  to  separate 
the  scattered  field  into  up-  and  down-going  P.  SV  and  S1I  waves  so  that  rough  interface 
scattering  coefficients  can  Ire  calculated.  To  test  the  perturbation  approximat  ion.  we  compare 
these  analytical  results  with  finite  difference  solutions  generated  for  the  same  rough  interface 
models.  Finally,  we  discuss  the  features  of  the  three-dimensional  scattered  field. 

THREE-DIMENSIONAL  SCATTERING  FORMULATION 

A  fine-scale  blow-up  of  a  three-dimensional  rough  interface  is  shoxvn  in  Figure  1 .  The  irregular 
interface  is  described  by  z  =  h(.v,y).  and  has  a  downward  normal  U_(x,y).  It  separates  an 
upper  medium,  described  by  P  and  S  wave  velocities  r>i  and  0\  and  density  /q,  from  a 
lower  medium,  described  by  02,  02,  and  Pi-  The  essence  of  the  formulation  is  to  project 
displacements  and  stresses  on  the  two  sides  of  the  rough  interface  onto  a  planar  surface 
whose  depth  equals  the  mean  depth  of  the  rough  interface.  These  projected  fields  are  then 
expanded  in  a  perturbation  series  in  h  about  a  background  field  consisting  of  the  known 
planar  interface  solution.  This  procedure  result s  in  a  formulation  in  which  the  rough  interface 
scattering  problem  is  replaced  by  a  planar  interface  scattering  problem,  with  sources  along 
the  planar  interface  generating  the  rough  interface  component  of  the  scattered  field. 

The  general  form  of  the  elastic  wave  equation  which  is  valid  for  small  displacements  in 
the  absence  of  body  forces  is 


-  pu>2v,  -  Tjij,  i,j  e  y,  ~}  ( 1) 

where  a,  is  the  fth  component  of  the  displacement  vector,  rJt  is  the  i.j  component  of  the 
Cauchy  stress  tensor,  the  comma  denotes  differentiation  in  the  y-th  coordinate  direction,  w 
is  the  temporal  frequency,  and  the  Finstein  summation  convention  applies.  Throughout  this 
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paper,  the  Fourier  transform  in  x,  y ,  and  t  uses  an  implied  phase  factor  of  exp [ikx.r  +  ikyy  — 
iut).  For  an  isotropic  solid,  stresses  are  related  to  displacements  by  the  constitutive  relation 


T{j  —  Xufc  k&ij  T  “h  (  —  ) 

where  Si:  is  the  Kroniker  symbol  and  A  and  p  are  the  Lame  parameters.  Equations  ( 1  j  and 
(2)  can  be  expressed  as  a  6  x  6  matrix  wave  equation  in  the  form 
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with  (  =  4/x(A  +  fi)/( A  +  2p). 

Welded  boundary  conditions  at  the  rough  interface  require  continuity  of  displacement 
and  traction  at  each  point  on  the  interface.  These  tractions  must  be  measured  with  respect 
to  a  plane  tangent  to  the  interface,  and  are  given  by  Tj  =  cr]knk-  The  downward  unit  normal 
nk  is  defined  by 
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A  new  vector  quantity  which  is  continuous  at  the  rough  interface  is  defined  to  be 
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Writing  (6)  in  a  form  which  explicitly  expresses  the  h  dependence,  this  transformation 
takes  the  form 

r(h)  =  (L  +  h,xgx  +  hygy)r(h)  (7) 

where  r  and  r  are  evaluated  along  the  rough  interface,  /  is  the  identity  matrix  and  Q  T  and 
Q_  y  are 
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The  interface  boundary  conditions  may  then  be  expressed  in  the  form 

r^(h)  =  f  <‘)(/j)  (10) 

where  superscripts  indicate  the  respective  media. 

Recalling  that  the  goal  is  to  relate  the  scattered  field  to  the  background  field  (scattered 
from  the  mean  planar  surface),  r  (h)  is  extrapolated  to  the  mean  planar  surface  by  the  power 
series  expansion 

r(h)  =  r  (0)  +  hr_\3z\ 0)  +  ,«(0)  +  . . .  (ID 

Making  use  of  the  wave  equation  (3),  (11)  is  reduced  to  the  form 

r(h)  =  (l  +  hA  +  f!~A2  +  ...)r(0).  (12) 

So  far  the  formulation  is  exact  as  long  as  the  series  in  (12)  converges.  It  is  easy  to 
demonstrate  convergence  for  the  case  of  a  planar  interface.  This  series  then  converges  to 

r(h)  =  eh=r  (0)  .  (13) 

The  meaning  of  eh=  is  easily  expressed  in  the  (fcx,  ky)  Fourier  transform  domain.  For  constant 
h ,  the  transform  is  done  by  replacing  dx  and  dy  in  A  by  ikx  and  iky ,  respectively,  and  formally 
summing  the  series  in  (12).  The  result  is  the  standard  form  for  the  propagator  matrix  (Aki 
and  Richards,  1980,  p.  275),  which  is  an  exact  extrapolation  operator  and  the  basis  of  the 
propagator  matrix  method  for  formulating  wave  equation  solutions  in  plane  layered  media 
(Kennett,  1983).  It  is  clear  that  a  simple  form  for  (12)  cannot  be  obtained  when  h  is  allowed 
to  vary. 

After  the  displacement-stress  vectors  along  the  rough  interface  have  been  extrapolated 
to  the  mean  planar  interface,  they  are  expanded  in  a  perturbation  series  about  ro(0).  the 
displacement-stress  vector  at  the  interface  of  the  background  planar  interface  model: 
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(M) 


r^>(0)  =  r  o(0)  4-  hr  (,j)(0)  +  h2r}j\ 0)  +  . . . 

Since  rVO)  is  the  displacement-stress  vector  for  the  background  planar  interface  model,  it  is 
continuous  across  the  boundary  and  needs  no  superscript.  The  higher  order  terms  reflect  the 
influence  of  the  rough  interface.  Combining  (7),  (12),  and  (14),  each  side  of  the  boundary 
condition  expressed  in  (10)  is  written  as  (omitting  the  superscripts) 

r_(h)  =  (£  +  4-  h  yQ_.  ) 

h2 

•  (X  +  hA  +  2?42  +  •••)  (b'>) 

( £  o  +  hi_  i  +  /i2rj  -4  .  .  .) 

~  (L  4  h  ,zQ_z  4  h'VQ  y)r0  -f  hA  r  g  4-  hr_  \  (Id) 

4  0(h?)  4-  0{hh<x)  4-  0(hhxy) 

where  O(-)  denotes  order  of  accuracy. 

Applying  the  boundary  condition  (TO)  to  (14)  results  in 


«ri2,-d")  =  h(Aw  -4'2,fco  +  /..,(£t',-£t2’)!;o  (IT) 

+  * .»(£', ”  -  £™)eo  - 

The  right-hand  side  of  (17)  would  be  zero  in  a  planar  interface  model.  Discontinuities  in  the 
displacement-stress  vector,  such  as  occurs  here,  represent  sources(Aki  and  Richards,  1980, 
p.  38).  Thus,  the  effect  of  a  rough  interface  is  to  appear  to  add  sources  along  the  mean 
planar  interface.  These  sources  will  be  designated  by  s  =  h{r  ^  —  r  j1').  This  use  of  sources 
to  represent  material  deviations  from  a  background  model  is  shared  with  standard  Born 
theoretical  developments.  The  mapping  of  heterogeneity  into  source  terms  is  also  used  in 
exact  formulations  based  on  Huygen’s  principle  (Paul  and  Campillo,  1988). 

Fourier  transforming  (17)  to  the  kx,  ky  domain, 
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whom  L_  is  defined  by 
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Definition  of  the  Scattering  Kernel 

The  perturbation  method  developed  above  will  now  be  used  to  examine  the  dominant  features 
of  the  three-dimensional  field  scattered  from  a  rough  interface.  For  plane  wave  sources  it 
is  possible  to  define  a  factorization  of  the  scattered  field  into  a  product  of  the  wavenumber 
spectrum  of  the  interface  and  a  function  that  is  called  the  scattering  kernel.  This  scattering 
kernel  is  independent  of  the  interface  roughness  function.  For  a  plane  wave  source  of  the 
form 


Lo(kx,  k'y)  =  4n2rp6{kx  -  kpx ,  k'y  -  kpy ), 


equation  (18)  reduces  to 


s(kx,ky)  =  h(kx  ~  kp,ky  -  kp)L(kx,ky-kp,kp)r0(kp,kp] 
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I  he  scattered  field  source  term  in  this  case  is  separated  into  a  part  associated  with  a  pniiii 
ular  interface  roughness,  h ( kT  —  A:£,  A-y  —  A-£),  and  a  part  associated  wit h  the  general  scattering 
properties  of  the  medium.  L.[kT,  ky;  k?,  A-j))r  0(  k?,  A-£).  The  latter  part  is  designated  as  tic 
scattering  kernel.  Knowledge  of  the  scattering  kernel  allows  one  to  evaluate  the  scattering 
potential  of  a  material  contrast  independent  of  any  particular  interface. 

The  transmission  scattering  kernels  for  the  two-dimensional  mode!  in  Figure  ts  are  given 
in  figure  2.  Superimposed  on  these  plots  is  the  Fourier  transform  ■>!  the  interface,  lie 
scattered  field  is  simply  the  product  of  the  two  curves.  From  these  plots  it  is  clear  that 
a  smaller  correlation  length,  and  hence  a  wider  Gaussian  in  the  transform  domain,  would 
result  in  large  amplitude  cusps  for  large  scattering  angles  in  P  and  S.  For  the  transmitted  P 
wave,  t  he  scat  tered  wave  amplitude  increases  for  scat  ter  mg  angles  larger  than  the  P  reflect  ion 
critical  angle.  This  amplitude'  boost  at  large  angles  may  be  seen  in  the  P  and  S  reflection 
coefficients  for  models  D  and  K  in  Figures  13(d  and  e).  This  effect  has  also  been  demonstrated 
by  Levander  and  Hill  (1985)  using  finite  difference  methods  and  by  by  Paul  and  Campillo 
( 19*8)  using  boundary  integral  equation  methods. 

Describing  Reflection  and  Transmission 

1  he  source  term  given  by  (18)  generates  P  and  S  waves  above  and  below  the  interface. 
To  determine  the  displacement  coefficients  of  these  waves  requires  a  relation  between  the 
displacement-stress  vector  and  the  up-  and  down-going  P,  SV.  and  SH  wave  components, 
the  form  of  which  is  given  by 

r  =  F_b  .  (22) 

where  b  -  [PS  T  P  S  T]1 ,  •  and  •  denote  down-  and  up-going  waves,  respectively,  and  1\ 
S.  and  T  are  the  displacement  coefficients  of  P,  SV',  and  SH  waves,  respectively.  Using  the 
reflection  coefficient  sign  conventions  of  Aki  and  Richards  (1980).  £  is  given  by 
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with  K  —  yjk%  -f  k2,  7  =  yju2 /a2  —  K2,  and  v  =  \Ju;2  /  fl2  —  K2.  To  recover  scattered  field 
displacements  from  the  source  term  note  that 

4  =  r  (2)  -  r'l>  =  £<2>6(2)  -  £('>6(1)  (24) 

and  that  s  generates  no  down-going  waves  in  the  upper  medium  and  no  up-going  waves  in 
the  lower  medium.  Hence, 
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The  wave  displacement  coefficients  generated  by  s  are  given  by 


b=t\  ('-Mi ) 

The  inverse  of  exists  for  all  values  cf  {kx.ky)  except  ( kT ,  ky )  =  (0,0).  At  this  point  in 
the  k.  plane  the  scattered  waves  are  propagating  normal  to  the  mean  planar  surface,  and  S\ 
waves  are  indistinguishable  from  SH  waves.  For  example,  consider  an  S  plane  wave  traveling 
in  the  c  direction  with  particle  motion  in  the  x  direction.  If  the  wave  direction  is  slightly 
perturbed  in  the  x  direction  it  becomes  an  SV  wave.  A  perturbation  in  the  y  direction  makes 
it  an  SH  wave.  These  distinctions  are  true  even  when  the  perturbations  are  infinitesimal. 
Hence,  to  remove  the  singularity  of  F,  an  arbitrary  naming  convention  must  be  adopted  for 
vertically  propagating  S  waves.  In  this  paper,  the  x  component  of  such  waves  is  labeled  SV. 
and  the  y  component  is  labeled  SH.  The  F_  matrix  for  vertical  propagation  is 
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COMPARISON  WITH  FINITE  DIFFERENCE 

The  simple  form  of  the  scattered  field  source  term  (18)  was  made  possible  by  the  truncation 
of  (15)  to  yield  a  low  order  perturbation  approximation  (16).  The  error  resulting  from  the 
exclusion  of  higher  order  terms  is  difficult  to  evaluate  analytically.  It  is  possible,  though,  to 
determine  bounds  on  the  domain  of  validity  of  this  approximation.  Kennett  (1972)  derived 
two  conditions  on  the  model  which  must  hold  in  order  for(18)  to  be  valid.  The  first  condition 
constrains  the  scattered  field  to  be  much  weaker  than  the  background  field,  a  requirement 
for  the  single  scattering  approximation  to  apply.  This  condition  is  expressed  by  the  relation 

■l(kj:,ky-,h  -  0)  <  max  |r0(A-r,  ky)\,  (28) 

fcj.fcy 

where  s  is  the  scattered  field  source  term  defined  by  (18).  Kennett  (1972)  reduced  this  to  a 
simpler,  but  stricter,  form  by  replacing  the  convolution  integral  in  (IS)  by  an  upper-hound 
approximation,  yielding 


KqujL 

— - — 7/12  max  |/i(x)j  <C  1 ,  max  \h  X(x)|  <C  1 , 
7rpi  1  T 


where 


r/12  =  max 


|Qipi  —  a2p  2 


I  P\  +  t>2/>2 


P\P\  ~  ihpi 
th  Pi  +  $2p2 


L  is  the  periodicity  length  of  the  rough  interface,  and  the  wavenumber  spectrum  of  the 
incident  field  is  bounded  by  |A;|  <  /c0.  The  second  condition  is  that  the  background  field  must 
not  contain  wavenumbers  so  close  to  grazing  incidence  that  shadow  zones  form.  Shadow  zones 
will  be  avoided  if  the  radius  of  curvature  of  the  interface  is  much  longer  than  a  wavelength. 
Such  waves  will  then  propagate  as  guided  waves  along  the  interface.  This  condition  is 
expressed  by 


— - -  max  \h  TT\  <C  1,  (30) 

+  ‘  ' 

where  kz  is  the  vertical  wavenumber  component  associated  with  the  maximum  horizontal 
component  /c0. 

These  conditions  are  not  satisfactory  for  practical  use,  however.  Approximations  used  in 
the  derivation  of  (29)  result  in  a  much  stricter  bound  than  is  necessary,  and  the  practical 
limit  imposed  by (30)  needs  to  be  better  defined.  In  order  to  empirically  construct  more 
realistic  constraints  on  interface  roughness,  reflection  and  transmission  coefficients  obtained 
using  the  perturbation  method  described  above  will  be  compared  to  those  derived  from  two- 
dimensional  finite  difference  solutions.  By  comparing  results  for  a  range  of  interface  height 
and  slope  statistics,  using  a  normal  incidence  plane  wave  source,  the  domain  of  validity  of 
the  perturbation  method  can  be  explored.  Three-dimensional  model  comparisons  using  the 
finite  difference  method  were  not  possible  with  the  computational  resources  available  for  this 
study.  The  similarity  of  the  two-  and  three-dimensional  perturbation  formulations  allows  the 
results  of  the  two-dimensional  validity  check  to  be  applied  to  three-dimensional  models. 

The  finite  difference  algorithm  used  here  is  a  two-dimensional  second-order  formulation. 
To  summarize,  the  wave  equations  (1)  and  (2)  are  solved  in  the  time-space  domain  bv  replac¬ 
ing  the  time  and  space  derivatives  by  their  second-order  centered-difference  approximations. 
The  accuracy  is  improved  by  using  a  staggered  mesh  formulation  in  which  horizontal  and 
vertical  displacements  are  represented  on  separate  grids,  each  shifted  by  half  of  the  grid 
point  spacing  in  both  coordinate  directions  with  respect  to  the  other.  This  has  the  added 
benefit  of  increasing  the  grid-point  density  by  a  factor  of  two.  The  formulation  differs  from 
that  of  Virieux  (1986)  in  that  we  use  the  second-order  displacement /stress  wave  equations 
instead  o  the  first-order  velocity/stress  equations.  This  modification  improves  efficiency, 
since  the  final  solution  was  desired  in  terms  of  displacement.  The  stability  conditions  and 
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grid  dispersion  relations  are  identical  with  those  of  Yirienx's  formulation.  A  low-order, 
staggered- mesh  scheme  was  chosen  because  it  is  the  most  efficient  method  when  dense  grid 
point  spacings  are  necessary,  as  is  the  case  in  our  models  where  small  interface  irregularities 
must  be  represented.  Higher  order  schemes,  such  as  those  which  use  fourth-order  (Bayliss. 
et  al.)  or  Fourier  spatial  derivative  operators  ( Kosloff  cl  al..  1981;  Fornherg,  198s).  are 
generally  thought  to  be  more  efficient  than  second-order  schemes,  but  this  is  not  necessarily 
true.  The  dense  grid  point  spacing  used  in  modeling  small  interface  perturbations  puts  the 
second-order  derivative  operator  well  within  its  domain  of  acceptable  accuracv.  and  the  short 
operator  length  makes  it  t fie  most  efficient  scheme. 

The  finite  difference  method  is  known  to  be  accurate  when  the  wave  fields  and  the 
model  are  both  well  discretized.  Hence,  when  models  with  large  interface  irregularities 
are  compared,  the  finite  difference  method  will  be  used  as  the  standard.  On  the  otlu-i 
hand,  the  perturbation  method  is  accurate  for  models  with  small  interface  height  and  slope 
irregularities  and  will  he  used  as  the  standard  for  such  models.  1  he  comparison  wit1,  finite 
difference  solutions  will  proceed  in  two  parts.  First  it  is  necessary  to  show  that  the  range  of 
interface  irregularities  over  which  the  finite  difference  method  is  accurate  extends  into  l  Ik' 
range  over  which  the  perturbation  method  is  accurate.  This  will  he  done  by  showing  that 
the  finite  difference  method  results  match  the  perturbation  method  results  for  a  model  with 
very  small  height  and  slope  irregularities.  If  the  finite  difference  method  works  well  in  this 
case,  results  for  larger  interface  irregularities  will  also  he  valid  because  they  will  be  more 
accurately  represented  on  the  finite  difference  grid.  Finite  difference  solutions  will  then  he 
used  a  basis  for  comparison  with  perturbation  solutions  for  a  series  of  models  with  larger 
height  and  slope  irregularities  in  order  to  probe  the  limits  of  validity  of  the  perturbation 
approximation. 

Accuracy  of  the  Finite  Difference  Method 

All  results  for  scattering  from  a  rough  interface  in  this  paper  are  expressed  as  reflection 
and  transmission  coefficients.  The  procedure  for  deriving  these  coefficients  from  the  finite 
difference  solution  is  similar  to  that  described  in  Section  for  the  perturbation  solution.  In 
short,  equation  (22)  is  used  to  convert  r  into  b,  and  then  P  and  S  wave  amplitudes  in  b 
are  scaled  by  the  source  wave  amplitude.  The  displacements  and  stresses  in  r  are  computed 
by  the  finite  difference  method  in  the  time-space  domain  along  a  horizontal  linear  array  of 
uniformly  spaced  receivers  that  do  not  intersect  the  interface  at  any  point.  The  receiver 
array  should  be  between  the  source  and  the  interface  in  order  to  separate  incident  waves 
from  reflected  waves,  r  is  then  Fourier  transformed  into  the  uj-k  domain  1o  yield  a  form 
suitable  for  use  in  (22).  The  distance  of  the  receiver  array  from  the  interface  is  accounted  for 
by  the  z  term  in  the  definition  of  F_  in  equation  (23).  If  the  finite  difference  model  is  such 
that  reflections  from  the  edge  of  the  grid  will  arrive  within  the  seismogram  time  window. 
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absorbing  boundaries  must  be  implemented  with  attenuate  these  reflections  to  a  level  much 
smaller  than  the  reflection  and  transmission  coefficients  to  be  measured.  File  length  of  the 
array,  L ,  controls  the  wavenumber  sampling  increment,  AA'  =  2n / L,  and  should  be  set  to  the 
horizontal  dimension  of  the  finite  difference  grid  to  avoid  aliasing.  From  here  forward  L  is 
taken  to  be  equal  to  the  horizontal  dimension  of  the  grid.  For  a  point  source,  L  controls  the 
range  of  incidence  and  scattering  angles  that  are  within  the  receiver  array,  and  the  horizontal 
grid  size  should  be  large  enough  to  capture  the  incidence  and  scattering  angles  of  interest. 
Since  the  range  of  incident  and  scattering  angles  detected  is  also  dependent  on  t  he  distance  of 
the  source  and  the  array  from  the  interface,  it  is  best  to  put  the  source  and  the  receiver  array 
as  close  to  the  interface  as  possible  in  order  to  minimize  the  grid  size  and  the  seismogram 
length.  The  spacing  of  the  receivers  controls  the  maximum  representable  wavenumber  in  r , 
and  should  be  set  equal  to  the  grid  spacing  to  avoid  aliasing. 

The  reflection  and  transmission  coefficients  obtained  from  finite  difference  modeling  are 
be  compared  with  analytical  coefficients  in  the  case  of  a  planar  interface  model.  The  model 
is  shown  in  Figure  3.  The  source  has  an  IS  Hz  Ricker  time  function  of  the  form 

R(t)  =  (\  (31) 

where  oj0  is  the  primary  angular  frequency  of  the  wavelet.  The  source  is  implemented  as  a 
body  force  term  of  the  form 


/  =  m 


— (x  ~  x,)e~(l  x’f!2 
-( 2  -  z,)z-(z-z’)2!2 


(32) 


representing  an  explosion  source  smoothed  by  a  Gaussian  with  standard  deviation  equal  to 
the  grid  point  spacing.  This  smoothing  is  necessary  to  make  the  explosion  dipoles,  which 
contain  energy  at  infinite  wavenumbers,  representable  on  the  finite  difference  grid.  The  source 
isotropically  radiates  pure  P  waves.  It  is  located  20  grid  points  above  the  interface,  and  the 
two  arrays  of  receivers  are  located  10  grid  points  above  and  10  grid  points  below  the  interface. 
Reflection  and  transmission  coefficients  derived  from  the  finite  difference  seismograms  for 
this  model  are  shown  in  Figure  4  along  with  analytical  coefficients.  There  is  generally  good 
agreement  for  both  pre-  and  post-critical  waves.  The  small  disagreement  present  at  the 
larger  scattering  angles  results  from  the  finite  aperture  of  the  receiver  array. 

Another  check  on  the  accuracy  of  the  reflection  and  transmission  coefficient  results  is 
to  plot  the  coefficients  for  the  non-physical  waves  in  the  solution:  the  down-going  S  in  ihe 
upper  medium  and  the  up-going  S  and  P  in  the  lower  medium.  These  coefficients  are  shown 
in  Figure  5.  For  each  wavenumber,  the  coefficients  are  normalized  by  the  amplitude  of  the 
down-going  P  wave  source.  Hence,  the  down-going  P  wave  coefficient  has  unit  amplitude 
for  all  angles.  The  remaining  three  coefficients  plotted  here  arc  indicative  of  error  in  the 
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conversion  of  the  measured  finite  difference  displacements  and  stresses  into  wave  coefficients. 
In  general,  the  .mplitude  of  the  non-physical  waves  increases  with  increasing  scattering 
angle,  since  the  wavenumbers  corresponding  to  these  larger  angles  are  less  well  represented 
by  the  receiver  array. 

The  effect  of  changing  the  array  aperture  is  illustrated  by  the  results  of  another  finite 
difference  model  w’hich  is  identical  to  the  previous  model  (Figure  3}  except  that  the  grid 
dimensions  and  receiver  array  sizes  are  doubled.  The  spacing  between  receivers  in  each 
array  is  unchanged.  A  comparison  of  reflection  and  transmission  coefficients  derived  from 
this  model  with  analytical  results  is  shown  in  Figure  6.  The  non-physical  coefficients  for  this 
model  are  shown  in  Figure  7.  The  ringing  observed  in  Figure  4  is  amplified.  In  addition, 
the  seismogram  time  window  in  this  model  was  inadvertently  set  large  enough  to  allow  the 
upward-traveling  source  wave  to  reach  the  top  of  the  grid  and  wrap  around  to  the  bottom 
of  the  grid  (we  are  using  periodic  boundary  conditions)  to  intersect  the  lower  receiver  array. 
This  results  in  the  large  lower  medium  up-going  P  wave  coefficient  shown  in  Figure  7  at  small 
scattering  angles,  and  to  a  lesser  extent  affects  the  accuracy  of  the  lower  medium  down-going 
P  wave  coefficient  shown  in  Figure  6. 

The  finite  difference  and  perturbation  methods  will  now  be  compared  for  a  model  with 
very  small  height  and  slope  perturbations  shown  in  Figure  8.  The  source  is  a  downward 
propagating,  planar  P  wave  with  an  18  Hz  Ricker  time  function.  The  interface  roughness 
has  a  Gaussian  autocorrelation  function  with  a  correlation  length  of  L  =  100  m  and  an  rms 
height  deviation  of  a  —  16.7  m.  The  interface  is  periodic  with  a  period  of  2.70  km,  the  width 
of  the  model.  The  F ourier  transform  of  the  interface  function  which  has  been  discretized  for 
the  finite  difference  grid  is  shown  in  Figure  9a  as  a  function  of  horizontal  slowness  evaluated 
at  18  Hz.  The  value  is  set  to  zero  at  the  origin  in  order  to  have  zero  mean  interface  deviation. 
It  would  be  a  perfect  Gaussian  were  it  not  for  the  finite  difference  discretization.  Histograms 
of  interface  height  and  slope  are  shown  in  Figures  9b  and  9c.  The  largest  deviation  from  the 
mean  planar  surface  is  12  m,  or  15  percent  of  the  S  wavelength  in  the  upper  medium.  The 
largest  slope  is  0.16,  while  the  mean  and  standard  deviation  of  the  slope  are  zero  and  0.059, 
respectively.  Since  the  maximum  interface  height  and  slope  are  fairly  well  approximated  by 
twice  their  standard  deviations,  histograms  will  not  be  provided  for  the  remaining  rough 
interfaces  used  in  this  study.  A  comparison  of  the  reflection  and  transmission  coefficients 
derived  from  the  finite  difference  and  perturbation  methods  is  shown  in  Figure  10.  Agreement 
is  very  good,  with  both  amplitude  and  shape  predicted  well  by  the  perturbation  method. 
'‘Ringing”  appears  to  some  degree  on  all  of  the  finite  difference  coefficient  plots,  and  is  most 
appaient  on  the  S  wave  coefficient  plots.  Plots  of  the  non-physical  coefficients  for  the  finite 
difference  data  are  shown  in  Figure  11.  The  down-going  P  wave  in  the  upper  medium  has 
unit  amplitude  at  zero  angle.  It  is  clear  from  this  figure  that  the  maximum  error  is  less 
than  0.1  percent,  and  that  this  error  is  associated  with  the  S  wave  in  the  upper  medium. 


This  error  is  enough  to  explain  the  difference  between  the  finite  difference  and  perturbation 
results  for  the  reflected  S  wave  coefficient  at  large  scattering  angles. 


Domain  of  Validity  for  the  Perturbation  Method 

The  domain  of  validity  of  the  perturbation  method  will  be  explored  by  comparing  reflection 

and  transmission  coefficients  generated  using  the  perturbation  method  with  those  derived 

from  finite  difference  modeling.  Six  rough  interface  models  are  used  in  this  comparison,  each 

having  a  Gaussian  autocorrelation  function  and  a  uniform  random  phase.  1  he  interface 

functions  are  shown  in  Figure  12(a-f)  with  two  times  vertical  exaggeration.  RMS  interface 

slope  ranges  from  0.037  to  0.99  and  RMS  interface  height,  is  0.01  km  for  models  A- K  and 

0.015  Km  for  model  F.  Since  the  accuracy  of  the  perturbation  method  is  sensitive  to  the 

smallest  elastic  wavelength  in  the  model,  interface  height  is  measured  here  in  terms  of  S 

wavelengths  in  the  upper  medium  (Sj),  and  thus  the  relative  height  of  the  irregularities 

changes  as  the  frequency  changes.  The  bandwidth  of  the  18  Hz  Ricker  wavelet  time  function 

used  in  the  finite  difference  calculation  allowed  reflection  and  transmission  coefficients  to 

be  computed  at  three  frequencies  (9.93,  16.5,  and  26.5  Hz)  in  each  of  the  six  models,  for 

a  total  of  18  comparisons  for  each  coefficient.  In  these  comparisons,  RMS  interface  height 

ranges  from  0.069  to  0.28  Sj  wavelengths.  RMS  height,  correlation  length,  and  RMS  slope 

for  the  models  used  are  listed  in  Table  1.  The  material  parameters  are  Q]  =  2.50  km/s, 

/?i  =  1.44  km/s,  pi  =  1.00  g/cm3,  a2  =  3.00  km/s,  d2  =  1-73  km/s,  and  p2  =  1.00  g/cm3. 

The  source  is  a  normally  incident  planar  P  wave  in  layer  one.  These  comparisons  are  shown 

for  models  A-F  in  Figures  13(a-f).  The  reflection  and  transmission  coefficients  are  the 

/  /  \  \  \ 

displacement  coefficients  Px,  Si,  P2,  and  S2  are  normalized  so  that  Px  =  1. 

In  order  to  facilitate  the  comparison  of  the  perturbation  and  finite  difference  coefficients, 
the  curves  are  compared  by  finding  the  L2  norm  difference  defined  by 


||Cp<  —  C/j II 2  —  100  x 


ff[PM)  -  P/d(0.)]2M. 

f}  Pfd(0,)3d9, 


(33) 


where  0,  is  the  scattering  angle  and  the  P  coefficient  is  chosen  for  illustration.  The  differences 
between  method  results  are  plotted  for  each  of  the  four  coefficients  in  Figures  !4(a-d)  for 
constant  RMS  interface  height,  and  in  Figures  14(e-h)  for  constant  RMS  interface  slope,  f  he 
constant  height  plots  show  that  the  accuracy  of  the  perturbation  solution  is  not  significantly 
degraded  as  RMS  interface  slope  varies  within  the  range  tested.  This  L2  norm  result  is 
verified  in  the  coefficient  comparison  plots  in  Figures  13(a-f).  The  constant  slope  plots  show 
that  error  increases  rapidly  with  increasing  RMS  interface  height,  and  in  this  particular 
example  is  acceptable  for  values  of  RMS  height  less  than  about  0.1  Si  wavelengths. 


FEATURES  OF  THE  3-D  SCATTERED  FIELD 


Three-dimensional  scattering  examples  will  now  be  presented  for  the  model  shown  in  Fig¬ 
ure  15.  The  material  and  interface  parameters  are  identical  to  those  of  model  D  (see  Ta¬ 
ble  1  and  Figure  12)  with  the  Gaussian  auto-correlation  function  having  identical  correlation 
lengths  in  the  x  and  y  directions.  In  general,  a  Gaussian  rough  interface  whose  major  and 
minor  axes  coincide  with  the  x  and  y  axes  has  a  Gaussian  auto-correlation  function  of  the 
form 


ky)hm(kT,  ky) 


If  f.j  ■  V,  ty  2 

L~  * 


(•'Til 


where  Lx  and  Ly  are  the  major-  and  minor-axis  correlation  lengths.  Equation  (3-1)  can  be 
generalized  to  allow  roughness  trends  at  an  angle  0  to  the  x  axis  by  applying  a  rotation 
transformation  to  get 


h(kx,  ky)hm(kx ,  ky) 


j[L2 (fc*  coi(6)+ky  8in(0))2  +  L2  (-/;x  sir>(5)-t-fcu  cos(0))2] 


(35) 


In  the  first  example,  the  source  is  an  18  Hz  planar  SV  wave  propagating  downward  along 
the  z  axis  with  particle  motion  in  the  x  direction.  The  SV  wave  scatters  into  reflected 
and  transmitted  P,  SV,  and  SH  waves.  The  three-dimensional  P,  SV,  and  SH  transmission 
coefficients  for  this  model,  found  using  the  perturbation  method,  are  given  in  Figure  16. 
The  x  and  y  scattering  angles  in  these  figures,  referred  to  as  4>x  and  <j>y  in  the  text,  are 
measured  from  the  downward  2  axis  in  the  x-z  and  y-z  planes,  respectively.  They  are 
defined  by  4>x  =  sin-1  ( kxv/u> )  and  <t>y  =  sin-1  (kyv/u)  where  v  is  the  body  wave  velocity 
of  the  transmitted  wave  concerned.  Scattering  in  the  2-direction,  for  example,  is  given 
by  <t>x  —  0y  =  0.  The  transmission  coefficient  plots  cover  scattering  angles  in  the  range 
-90°  <  <f>x  <  90°  and  0°  <  (f>y  <  90°,  and  are  symmetric  about  the  <f>y  =  0  axis.  For  all 
three  scattered  wavetypes,  these  plots  show  that  scattering  is  maximal  in  the  direction  that 
conserves  source  particle  motion:  for  an  SV  source,  P  and  SV  are  maximally  scattered  in 
the  x-z  plane  (<j)y  =  0),  and  SH  is  maximally  scattered  in  the  y-z  plane  (<f>x  =  0).  This  effect 
is  exaggerated  by  the  use  of  the  single  scattering  approximation.  A  null  is  present  in  the 
plane  normal  to  maximal  plane;  the  y-z  plane  for  P  and  SV  waves  and  the  x-z  plane  for  SH 
waves.  An  alternative  display  of  the  three-dimensional  transmission  coefficients  is  a  cross- 
section  of  the  coefficients  for  all  scattering  angles  at  a  particular  azimuth,  where  azimuth  is 
measured  clockwise  in  the  x-y  plane  from  the  x-axis,  and  the  scattering  angle  is  defined  by 
0  =  sin-1  ( Kv/u> ).  Cross-sections  of  the  reflection  and  transmission  coefficients  for  the  same 
model  and  source  parameters  are  given  for  azimuthal  angles  of  10°,  45°,  and  80°in  Figure  17. 

Transmission  scattering  kernels  for  the  model  in  Figure  15  are  given  in  Figure  18.  Com¬ 
parison  with  the  SV  and  SH  transmission  coefficients  in  Figure  16  shows  that,  the  spatially 
band-limited  nature  of  this  particular  interface  damps  out  the  large  amplitude  features 
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present  at  large  scattering  angles.  The  P  transmission  coefficient  exhibit-  .<  hint  of  the 
cusps  present  in  the  kernel.  Consideration  of  another  interface  louehue--  fmu  » i« »s i  requires 
only  a  visual  superposition  of  the  new  interface  spectrum  with  the  heM.i lm  the  suin' 
SV  source  as  above  and  a  two-dimensional  rough  interface  with  vai  i.-r  the  r  direc¬ 

tion,  the  problem  is  truly  two-dimensional,  and  waves  are  scattered  ii.t.i  p  and  S\  If  the 
interface  is  rotated  90  degrees  so  that  variation  is  in  the  y  direction,  the  problem  I-  fully 
three-dimensional,  and  waves  are  scattered  into  P  and  Sll.  More  general  interfaces  whose 
auto-correlation  functions  are  described  by  (35),  or  perhaps  a  ton  Karman  or  exponential 
function,  are  handled  with  the  same  approach. 

The  second  example  is  the  same  as  the  first,  but  with  the  SV  source  replaced  by  a  P 
source.  The  three-dimensional  P  and  SV  transmission  coefficients  for  this  model  are  given  in 
Figure  19.  Since  the  source  particle  motion  and  the  Gaussian  auto-correlation  function  are 
azimuthally  invariant,  the  scattering  coefficients  are  also  azimuthally  invariant.  Deviations 
of  the  contours  from  circular  arcs  are  artifacts  of  the  contouring  program.  Reflection  and 
transmission  coefficient  cross-sections  for  this  model  are  given  in  Figure  20.  Note  that 
scattered  SH  waves  are  not  generated  in  this  example.  This  is  a  consequence  of  the  single¬ 
scattering  approximation.  For  a  normal  incidence  planar  P  wave,  a  minimum  of  two  bounces 
are  required  to  generate  SH. 

DISCUSSION  AND  CONCLUSIONS 

A  perturbation  method  has  been  presented  here  for  computing  three-dimensional  body  wave 
scattering  from  a  rough  interface.  Its  speed  and  simplicity  are  such  that  each  of  the  examples 
in  this  paper  was  generated  on  a  Macintosh  SE  computer  in  less  than  two  minutes.  Speed 
is  an  important  consideration  when  three-dimensional  modeling  is  necessary.  For  example, 
many  of  the  two-dimensional  finite  difference  computations  done  for  comparison  required  15 
Mbytes  of  core  and  23  hours  of  CPU  time  on  a  Vax  8800.  In  three  dimensions,  finite  difference 
solutions  for  these  models  would  be  beyond  our  resources.  For  the  class  of  irregular  interface 
models  with  small  RMS  height  deviations,  the  perturbation  method  is  a  useful  alternative. 

Since  interface  height  deviations  in  the  models  presented  here  are  small,  comparisons  of 
the  perturbation  method  with  the  finite  difference  method  were  preceded  by  careful  testing 
of  the  finite  difference  method  to  show  that  it  is  valid  for  the  small  deviations  used  in  these 
comparisons.  It  can  be  shown  that  in  the  limit  as  the  finite  difference  grid  sampling  interval 
goes  to  zero,  and  as  the  numerical  precision  of  the  computer  goes  to  infinity,  the  finite 
difference  solution  converges  to  the  exact  solution  as  0(A.r)  (Brown.  1984).  However,  in 
order  for  the  grid  sampling  interval  to  be  small  enough  that  the  finite  difference  solution  is 
sufficiently  close  to  the  exact  solution,  the  number  of  grid  points  describing  a  fixed  model 
must  increase  as  the  size  of  the  interface  height  perturbations  decreases.  Therefore,  the  core 
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size  and  speed  of  a  computer  are  constraints  in  the  minimum  interface  height  pert urhat ion 
that  can  be  accurately  modeled  using  finite  difference  method.  Of  course,  a  planar  boundary 
is  an  exception  to  this  limit.  Another  lower  limit  on  interface  perturbation  size  is  imposed  by 
numerical  precision.  As  the  size  of  interface  height  perturbations  decreases,  the  amplitude 
of  scattered  seismic  waves  decreases.  Since  these  scattered  waves  are  added  to  the  relatively 
large  planar  interface  response,  insignificant  numerical  noise  in  the  planar  interface  response* 
can  be  significant  when  compared  with  the  scattered  field.  This  is  probably  the  cause  of  the 
“ringing”  present  in  the  finite  difference-derived  scattering  coefficients.  Comparisons  of  finite 
difference-derived  scattering  coefficients  for  a  planar  model  against  an  analytical  solution, 
and  for  a  model  with  small  interface  height  and  slope  perturbations  against  a  perturbation 
method  solution,  show  that  this  “ringing”  tends  to  oscillate  about  the  exact  solution.  Hence, 
the  smoothed  finite  difference  scattering  coefficients  are  useful  for  comparisons. 

Numerous  two-dimensional  model  comparisons  of  finite  difference  and  perturbation  met  hod 
scattering  coefficients  are  given  in  Figures  13(a-f).  In  order  to  make  it  easier  to  spot  trends 
in  these  comparisons,  each  comparison  is  summarized  by  a  single  number,  the  b2  norm  of 
the  difference,  which  is  plotted  versus  RMS  slope  for  a  constant  RMS  height,  and  against 
RMS  height  for  a  constant  RMS  slope,  in  Figure  11.  There  are  two  major  trends  in  the 
L2  norm  results.  The  error  increases  strongly  with  increasing  RMS  height,  and  has  accept¬ 
able  levels  for  RMS  heights  of  less  than  about  0.1  Sj  wavelengths.  Also,  error  appears  to 
be  roughly  constant  for  increasing  RMS  slope  for  the  tested  range  of  from  0.037  to  0.09. 
Smaller  trends,  such  as  the  apparent  increase  in  accuracy  with  increasing  RMS  slope,  are 
misleading.  The  L2  norm  is  overly  sensitive  to  the  noise  in  the  finite  difference  solution  being 
used  as  the  standard  for  comparison.  This  is  apparent  when  the  scattering  coefficient  plots 
in  Figures  1 3  ( a  -  f )  are  examined  for  the  trend  seen  in  the  L2  norm  plots.  The  Lj  norm  was 
also  calculated  to  see  whether  it  is  less  sensitive  to  this  noise,  but  the  improvements  were 
minimal.  Ultimately,  all  trends  must  be  confirmed  by  the  scattering  coefficient  plots. 

The  perturbation  method  was  used  to  calculate  P,  SV,  and  SH  scattering  for  an  SV 
wave  normally  incident  on  a  three-dimensional  rough  interface.  It  was  shown  that  P  and  SV 
waves  scattering  in  the  azimuth  of  the  incident  SV  particle  motion  have  maximum  amplitude, 
and  that  no  P  or  SV  waves  are  scattered  in  the  azimuth  normal  to  the  maximal  direction. 
Scattered  SFI  waves  have  the  opposite  behavior,  scattering  maximally  in  the  azimuth  normal 
to  the  incident  SV  particle  motion.  A  scattering  kernel  was  defined  and  demonstrated  which 
allows  the  scattering  radiation  pattern  to  be  separated  into  a  part  controlled  by  the  material 
contrast  at  the  irderface  and  a  part  due  to  the  interface.  The  scattering  kernel  for  the  three- 
dimensional  scattering  example  shows  that  scattered  wave  amplitudes  tend  to  increase  in 
amplitude  for  scattering  angles  beyond  the  P  wave  critical  angle  for  the  background  planar 
interface  model. 

In  the  three-dimensional  scattering  example  it  was  shown  that  there  is  a  null  in  P  and 
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SV  scattering  in  the  y  direction,  and  a  null  in  SH  scattering  in  the  s  direction.  These  nulls 
do  not  exist  in  the  time-space  domain,  where  a  receiver  in  any  location  can  detect  waves 
traveling  in  all  directions.  Future  work  includes  an  extension  of  the  perturbation  results  to 
the  time-space  domain  using  the  discrete  wavenumber  method  (Bouchon,  1977)  in  order  to 
account  for  these  interference  effects. 
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Model 

Correlation 

Length 

■ 

RMS 

Slope 

RMS  Height 
in  Si  Wavelengths 
at  9.93  Hz 

RMS  Height 
in  Si  Wavelengths 
at  16.5  Hz 

RMS  Height 
in  Si  Wavelengths 
at  26.5  Hz 

A 

0.30 

0.037 

0.069 

0.11 

0.18 

B 

0.10 

0.10 

0.069 

0.11 

0.18 

C 

0.050 

0.20 

0.069 

0.11 

0.18 

D 

0.033 

0.30 

0.009 

0.11 

0.18 

E 

0.010 

0.99 

0.069 

0.11 

0.18 

F 

0.15 

0.10 

0.10 

0.17 

0.28 

Table  1:  Interface  parameters  for  the  rough  interfaces  shown  in  Figure  12. 
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Figure  1:  Ceometry  of  the  rough  interface  model.  The  interface  is  defined  by 
and  the  downward-pointing  unit  normal  to  this  surface  is  denoted  by  n_.  The 
is  defined  as  the  mean  planar  surface  for  the  rough  interface. 
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Figure  2:  Transmission  scattering  kernels  for  the  model  in  Figure  8.  The  Fourier  transform 
of  the  interface  is  superimposed. 
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Interface  Spectral  Amplitude 


Figure  3:  Planar  interface  model  used  to  compare  finite  difference  derived  reflection  and 
transmission  coefficients  with  analytical  solutions.  The  source  is  a  point  explosion  with 
an  18  Hz  Ricker  wavelet  time  function.  The  sampling  parameters,  Ax  =  A z  =  6.67  ru 
and  At  =  0.00156  s,  result  in  a  maximum  phase  dispersion  error  of  —1.1  percent  at  IS 
Hz.  The  source  is  located  20  grid  points  above  the  interface,  and  the  receiver  arrays  are 
located  10  grid  points  on  either  side  of  the  interface.  The  horizontal  receiver  interval  is 
two  grid  points.  All  waves  within  the  seismogram  time  window  are  contained  within  the 
finite  difference  grid,  eliminating  the  need  for  absorbing  boundaries. 
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Figure  4:  Comparison  between  P  and  S  wave  reflection  and  transmission  coefficients  gener¬ 
ated  by  finite  difference  and  analytic  methods  for  a  planar  interface.  The  model  is  shown 
in  Figure  3.  The  small  disagreement  present  at  the  larger  scattering  angles  is  ‘‘Gibb’s 
ringing”  that  results  from  the  finite  aperture  of  the  receiver  array. 
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Amplitude  Amplitude 


Down-going  P  Down-going  S 


Up-going  P  Up-going  S 


Figure  5:  These  plots  of  the  wave  amplitudes  incident  on  the  interface  provide  a  measure 
of  the  error  versus  angle.  For  each  wavenumber,  the  coefficients  are  normalized  by  the 
amplitude  of  the  down-going  P  wave  source.  Hence,  the  down-going  P  wave  coefficient 
has  unit  amplitude  for  all  angles.  The  remaining  three  coefficients  plotted  here  are  non¬ 
physical,  and  are  indicative  of  error  in  the  conversion  of  the  measured  finite  difTerrme 
displacements  and  stresses  into  wave  coefficients. 
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Transmission  Coefficient  Reflection  Coefficient 


Reflected  P  Reflected  S 
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rigure  6:  Comparison  between  P  and  S  wave  reflection  and  transmission  coefficients  gener¬ 
ated  by  finite  difference  and  analytic  methods.  The  finite  difference  model  is  identical 
to  the  previous  model  (in  Figure  3)  except  that  the  grid  dimensions  and  recever  array 
sizes  are  doubled.  The  receiver  spacing  is  unchanged.  The  only  effect  is  to  double  the 
aperture  of  the  receiver  array.  The  result  is  to  amplify  the  "Gibb’s  ringing”  effect  of  the 
Figure  4  results. 
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Figure  7:  Amplitudes  of  waves  incident  on  the  interface  corresponding  to  the  scatter* 
amplitudes  shown  in  Figure  6.  These  are  provided  as  a  measure  of  the  error  versu 


Figure  8:  Two-dimensional  rough  interface  model  used  to  compare  a  perturbation  solution 
with  a  finite  difference  solution.  The  source  is  a  normally  incident,  18  Hz,  plane  wave. 
The  rough  interface  has  a  Gaussian  autocorrelation  function  with  a  correlation  length  of 
L  —  100  m  and  an  rms  height  deviation  of  a  =  16.7  m.  The  interface  is  periodic  with  a 
period  of  2.70  km,  the  width  of  the  model  above.  The  sampling  parameters,  Ax  =  6.67 
m,  A z  —  2.22  m,and  A t  —  0.00156  s,  result  in  a  maximum  phase  dispersion  error  of  —1.1 
percent  at  18  Hz.  The  two  receiver  arrays  are  located  as  close  to  the  interface  as  is  possible 
without  intersection.  The  horizontal  receiver  interval  equals  Ax.  All  waves  within  the 
seismogram  time  window  are  contained  within  the  finite  difference  grid,  eliminating  the 
need  for  absorbing  boundaries. 
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Figure  9:  Properties  of  of  the  interface  shown  in  Figure  8:  (a)  Fourier  transform,  where  the 
horizontal  slowness  is  evaluated  for  18  Hz.  and  histograms  of  (b)  interface  height  and  (c) 
interface  slope. 
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Figure  10:  Comparison  between  P  and  S  wave  reflection  and  transmission  coefficients  goner 
ated  from  the  finite  difference  and  perturbation  methods  for  the  model  shown  in  Figure  s 
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Figure  11:  Amplitudes  of  waves  incident  on  the  interface  corresponding  to  the  scattered 
wave  amplitudes  shown  in  Figure  10.  These  are  provided  as  a  measure  of  the  error  versus 
angle.  The  down-going  P  wave  in  the  upper  medium  has  unit  amplitude  at  zero  angle. 
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Figure  12:  Interface  functions  used  in  comparison  of  reflection  and  transmission  coeffi¬ 
cients  derived  from  the  perturbation  method  with  those  derived  from  the  finite  difference 
method.  The  interfaces  have  a  Gaussian  autocorelation  function  with  RMS  height,  cor¬ 
relation  length,  and  RMS  slope  for  each  interface  listed  in  Table  1. 
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Figure  13a:  Comparison  of  reflection  and  transmission  coefficients  derived  from  tin*  |» ■!’ 
bation  and  finite  difference  methods  for  model  A.  The  parameters  of  the  rough  int.  ?■  ■  ■ 
are  given  in  Table  1  and  the  interface  function  is  illustrated  in  Figure  12. 
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Figure  13d:  Comparison  of  reflection  and  transmission  coefficients  derived  from  the  pertur 
bation  and  finite  difference  methods  for  model  D.  The  parameters  of  the  rough  inter fao 
are  given  in  Table  1  and  the  interface  function  is  illustrated  in  Figure  12. 
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Figure  14:  L2  norm  comparison  of  finite  difference  and  perturbation  method  derived  refle(  t  ion 
and  transmission  coefficients,  (a-d)  RMS  interface  height  is  a  constant  0.01  km  and 
RMS  slope  varies  from  0.037  to  0.09.  RMS  interface  height  for  three  frequencies  <  ,m 
he  expressed  as  0.009.  0.11,  and  0.18  S|  wavelengths,  (r  h)  RMS  interface  slope  is  ,, 
constant  0.1  and  with  RMS  interface  height  varies  from  0.069  to  0.28  Si  wave!r’>g<  hs. 
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Figure  15:  Three-dimensional  rough  interface  scattering  model.  The  auto-correlation  func¬ 
tion  of  the  interface  is  a  two-dimensional  Gaussian  with  x  and  y  correlation  lengths  of 
2.4  Si  wavelengths  (0.19  km),  an  RMS  height  of  0.125  Si  wavelengths  (0.010  km),  and 
an  RMS  slope  of  0.30.  Contours  and  axes  are  labeled  in  kilometers. 
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Figure  17:  Cross-sections  of  the  transmission  coefficients  in  Figure  16  along  10°,  45°,  and 
8Q°azimutiis.  Reflection  coefficients  for  the  same  model  are  also  shown. 
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Figure  20:  Cross-sections  of  the  transmission  coefficients  in  Figure  19  along  10°,  15 
S0°azimuths.  Reflection  coefficients  for  the  same  model  are  also  shown. 
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